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We have used an exact general relativistic model structure within a FRW cosmological background 
based on a LTB metric to study the gravitational lensing of a cosmological structure. The integration 
of the geodesic equations turned out to be a delicate task. We realized that the use of the rank 8(7) 
and 10(11) Runge-Kutta numerical method leads to a numerical effect and is therefore unreliable. 
The so-called semi-implicit Rosenbrock method, however, turned out to be a viable integration 
method for our problem. The deviation angle calculated by the integration of the geodesic equations 
for different density profiles of the model structure was then compared to those of the corresponding 
thin lens approximation. Using the familiar NFW density profile, it is shown that independent 
of the truncation details the thin lens approximation differ substantially from the exact relativistic 
calculation. The difference in the deflection angle for different impact parameters may be up to about 
30 percent. However, using the modified NFW density profile with a void before going over to the 
FRW background, as required by an exact general relativistic model, the thin lens approximation 
coincides almost exactly with the general relativistic calculation. 

PACS numbers: 98.80.Jk, 98.62.Js, 98.62.Ck , 95.35.+d 



The thin lens (Th-L) approximation in the gravi- 
tational lensing is the prevailing method to estimate 
cosmological parameters and the mass of large scale 
structures leading to dark matter and dark energy 
contents of the universe pQ, |2]. The current view is 
that this Th-L approximation is accurate enough at the 
cosmological scales where we are faced with very weak 
gravitational fields and potentials. There has already 
been attempts to compare the Th-L approximation 
with the integration of null geodesies in a perturbed 
cosmological background ([SHH], see also [5] and the 
references there). However, a full general relativistic 
calculation based on an exact model is still missing. 
There are two sources of misinterpretation of astro- 
physical phenomenon in a weak gravity environment, 
depending on the local or quasi-local phenomena under 
consideration. In the case of local phenomena the 
familiar perturbation theories maybe valid to some 
extend. There are already detailed studies on this 
subject (see [B], [7], [5]). However, if quasi- local phe- 
nomena or structures come into play we may encounter 
counter-intuitive effects not detected in the perturba- 
tional approach to the weak field limits. The definition 
of quasi-local mass in general relativity is one of these 
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issues which has been extensively studied in general 
relativity [5]. We have already shown numerically how 
different various quasi-local mass definitions of a general 
relativistic structure may be [10]. Another quasi-local 
effect relevant to the gravitational lensing is how a 
spherically symmetric structure is matched to a FRW 
background. Such a general relativistic matching is only 
possible through an underdensity region or a void [llj : 
a fact not realized in the post-Newtonian approaches 
or cosmological perturbations relevant to lensing, and 
missed in all studies comparing the Th-L approaches to 
a more exact general relativistic lensing calculation. 



We are interested in the exact general relativistic lens- 
ing by an exact solution of Einstein Equations repre- 
senting a cosmological structure defined by a spherically 
symmetric overdensity structure within a FRW universe. 
There is already an exact general relativistic model struc- 
ture within an FRW universe based on a Lemaitre, Tol- 
man and Bondi (LTB) metric |13t415j representing an in- 
homogeneous cosmological model with a structure at its 
center [TUj. Choosing such a model for an extended spher- 
ical lens, we study the gravitational lensing in a dynami- 
cal cosmological background in the framework of general 
relativity by integrating numerically the null geodesic 
equations to obtain the deflection angle. The result is 
then compared with the corresponding Th-L approxima- 
tion to understand the accuracy of this technology and 
its possible flaws in interpreting the structure and the 



mass of cluster of galaxies. The effect of the cosmolog- 
ical constant in the lensing is negligible in small scales 
we are considering [T2j and only effect the cosmological 
distances which we will take into account. That is why 
we have neglected the cosmological constant in our exact 
model to avoid unnecessary complexities 
Take a spherically symmetric cosmological structure in a 
FRW matter dominated universe with the density p(r 7 t) . 
This is modeled by a LTB solution of the Einstein equa- 
tions which is written in the comoving coordinates as 
(G = l,c=l) 



ds 2 



satisfying 



-dt 2 + X 2 (r, t)dr 2 + R 2 (t, r)dtt 2 . 
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Here M and E are integrating functions, where dot and 
prime denote partial derivatives with respect to the co- 
ordinates t and r respectively. Equation Q has three 
different analytic solution, depending on the value of E. 
The solution for negative E we are interested in is given 
by 



R 

rj — sin rj 



M , 

-^-(l-COST?), 



{-Ef / 2 
M 



<t-t b (r)). 



(5) 



The solution has three free functions: £f,(r), E(r), and 
M(r). Given that the metric is covariant under the 
rescaling r — > f(r) one of these functions may be fixed. 
The geodesic equations may be written in the arbitrary 
plane of 6 — | due to the spherical symmetry: 
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where A is an affine parameter. Equation ([8]) expresses 
the conservation of the angular momentum: 



L = R 2 ^= Const. 
dX 



(9) 



We are interested in the light-like geodesies. From the 
metric we obtain the light-like condition in the form 



= X 2 



R z 



dX 



(10) 



These partial non-linear differential equations can not be 
solved analytically. To integrate them numerically one 
has to specify the three functions M(r),ti,(r), and E{r) 
and all derivatives of the metric functions, using a proce- 
dure proposed in [T!)l[T?] . We start with a generic density 
profile and specify it at two different times t\, t2 as a func- 
tion of the coordinate r. Now, the numerical procedure is 
based on the choice of r-coordinate such that M(r) = r. 
This is due to the fact that M(r) is an increasing func- 
tion of r. Therefore, E and % become functions of M. 
For the initial time we choose the time of the last scat- 
tering surface: t\ ~ 3.77 x 10 5 yr. The initial density 
profile should show a small over-density near the center 
imitating otherwise a FRW universe. Therefore, we add 
a Gaussian peak to the FRW background density pt- We 
know already that having an over-density in an otherwise 
homogeneous universe needs a void to compensate for the 
extra mass within the over-density region. Therefore, to 
model this void we subtract a wider gaussian peak: 



p(R,h) = Pb(tx) 



5\e 



-h 
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where 5\ is the density contrast of the Gaussian peak, Rq 
is the width of the Gaussian peak, and R\ is the width of 
the negative Gaussian profile. The mass compensation 
condition leads to an equation for b\. For the final time 
we choose the time when our null geodesy has the near- 
est distance to the center of our model structure. For 
instance if we set our lens at the redshift z ~ 0.2 then 
t 2 ~ 6.98Gyr. 

The density profile we choose for the final time is the 
universal halo density profile (NFW) [TH] convolved with 
a negative Gaussian profile to compensate the mass plus 
the background density at that time: 



p(R, t 2 ) = (pnfw - b 2 Pb{h)) e 



where 



and 



Pnfw — Pfc(M 
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In our numerical calculation we will use typical NFW 
values R s — 0.5Mpc and c = 5 for a galaxy cluster. Note 
that at the time t 2 a black hole singularity covered by an 
apparent horizon has already been evolved. Therefore, 
the NFW profile has to be modified and a black hole mass 
greater than a minimum value has to be added to it at the 
center. This physical fact is reflected in a shell crossing 
singularity if we take the familiar NFW profile similar to 
that assumed for the time t\ . The mass we have assumed 
for this black hole singularity is about one thousandth of 





l.OE+11 l.OE+13 l.OE+15 l.OE+17 l.OE+19 l.OE+21 l.OE+23 l.OE+25 
M/M a 



1.0E-04 1.0E-03 l.OE-02 l.OE-01 l.OE+00 l.OE+01 l.OE+02 1.0E+03 1.0E+04 
R(Mpc) 



FIG. 1: E as a function of M for a cluster with NFW density 
profile. M is given in the unit of the Sun mass. 



FIG. 3: Density profile for a cluster. The dot line corresponds 
to the familiar NFW profile and the solid line corresponds to 
the modified NFW with a void. 
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The integration is done from the observer to the source 
at a specific redshift. Assuming there is no lens, the 
model reduces to a homogenous flat FRW universe and 
the geodesies are straight lines (in comoving coordinates) 
allowing us to determine the angle between the source 
and the lens (/3 in Fig. [4|: 
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FIG. 2: tb as a function of M for a cluster with NFW density 
profile, tb is given in the unit of 3.263Gj/r. 



the mass up to the R s and equal to 5.66 x 10 11 M Q . Figs. 
[T] and [2] shows the LTB functions E and tb as a result of 
these boundary assumptions. Using these LTB functions, 
the density profile of our model structure is obtained and 
depicted in Fig. [3] 

To solve these equations we have to specify four initial 
conditions taking into account the light-like condition 
(10). The freedom of choosing the affine parameter 
reduces the initial conditions to three. Now, the in- 



tegration of the geodesies happens by a backshooting 
procedure. Our initial conditions are taken to be the 
time of observation, distance of the observer to the lens 
expressed in terms of the redshift of the lens at the time 
of the observation, and angle between the line of sight 
to the image of source and the line of sight to the lens 
{6 in Fig. gj: 



where <f> e is the OLS angle, r is the comoving distance 
of the observer, and r e is the comoving distance of the 
source from the center of coordinate system in the ab- 
sence of lens at the time t e . From the geodesic equations 
the t e is given by 
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We then write the lens equation and determine the 
deflection angle 7: 



7 = {9 - p) 
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D L s 



(18) 



where we have assumed that the presence of the lens has 
not a significant effect on the distances and we may use 
the corresponding FRW ones. 

The validity of the numerical method chosen to 
integrate such complex system of partial differential 
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FIG. 4: GL diagram: O is observer, S is source, S' is image 
in source plane, L is lens and 7 is deflection angle. 



equations is a delicate issue. We first started with the 
familiar Runge-Kutta adaptive step size algorithm with 
proportional and integral feedback (PI control) JT5] in 
which the step size is adjusted to keep local error under 
a suitable threshold. We started with the so-called 
embedded Runge-Kutta of the rank 5(4). It turned out, 
however, that its accuracy is too low. Therefore, we 
tried the rank 8(7) and then the rank 11(10) algorithm. 
The difference between these two last ranks, however, 
turned out to be marginal and below one percent. Given 
the time-consuming rank 11(10) algorithm, we preferred 
to use the rank 8(7) one. Now, as a fist test for the 
accuracy of this numerical method we tried the trivial 
example of the LTB model, namely the FRW case, 
expecting a null result. The result was a non-negligible 
deflection angle of the order of few milliarcseconds. 
Suspecting to face a numerical effect, and trying to 
understand the numerical algorithm and the source of 
this numerical effect, we continued to calculate a more 
concrete and non-trivial LTB case. The result for the 
rank 8(7) Runge-Kutta numerical method applied to a 
structure with a compact density profile did agree with 
the thin lens approximation. However, in the case of a 
more diffuse density profile the result showed a deflection 
angle up to an order of magnitude higher than the thin 
lens approximation. We did interpret this result as a 
sign not to trust the Runge-Kutta method and turned 
to an alternative numerical method! 

The root of this numerical deficiency could be due to 
the term -3? in our equations, which is almost zero in 
the most part of the path of the light ray and changes 
suddenly to 7r in the vicinity of the lens. This is a 
well-known phenomenon in the numerical method of 
integrating differential equations called as "stiff" [20] , 
The characteristic property of such stiff equations is 
the presence of two quite different scales. In our case 
we have on one side the cosmological distance scale of 
the source relative to the lens and the observer, and on 
the other side the scale of the structure or the nearest 



distance of the ray to the lens. Realizing this stiffness 
property, we turned to the so-called semi-implicit 
Rosenbrock method of the numerical integration of 
partial differential equations [T^l [2D] . As a first test we 
calculated again the trivial case of a FRW model which 
gave an acceptable null result. We, therefore, decided to 
integrate our geodesic equations using the semi-implicit 
Rosenbrock method instead of the Runge-Kutta one. 

The null geodesies equations of our exact general rela- 
tivistic structure model is now integrated using the mod- 
ified NFW density profile with a void before matching to 
the background FRW universe to obtain the deflection 
angle. Note that the density in the NFW density pro- 
file is taken to be the oversdensity in an otherwise FRW 
model, namely p — pt- However, for the Th-L approxima- 
tion we have used two different density profiles namely 
the familiar one and the modified one with a void be- 
fore matching to the background density. In the case of 
familiar NFW density profile without a void, the corre- 
sponding equations can be integrated analytically to give 
the deviation angle (23] [25] : 
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Assuming the same modified NFW profile as in general 
relativistic case for the Th-L approximation we have also 
calculated the deflection angle applying the lens equation 
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where ^/ is the lens potential. 

The result for the three cases, the exact general rela- 
tivistic model with our modified NFW profile, thin lens 
approximation using the modified NFW with void, and 
the thin lens approximation using the familiar NFW 
without a void is depicted in Fig. [5j Obviously the two 
cases of the thin lens approximation with the modified 
NFW density profile including the void and the LTB ex- 
act method almost coincide. The thin lens approximation 
with the familiar density profile without a void, however, 
differ from the exact LTB model. The difference in the 
deviation angle can be more than 30 percent depending 
on the impact parameter. The difference between the 
exact general relativistic LTB model and the thin lens 
approximation is due to the absence of the void in the 
familiar NFW profile used in the literature. To see the 
implications of the NFW parameters in this difference 
we have also calculated the deviation angle for different 
NFW profiles, with and without void. The result is de- 
picted in the Fig. [19] We see again that the Th-L approx- 
imation using different modified NFW profiles including 
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FIG. 5: Deviation angle for three cases: the general relativis- 
tic result is indicated by plus points; the thin lens approxima- 
tion using our modified NFW is shown by the continuous line; 
and the dashed line is for the familiar NFW profile without 
the void (formula ( 19 1 ) . 
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a void almost coincide with the exact LTB model. Mod- 
els with the NFW profiles without void, however, differ 
substantially from the exact model. The difference is 
higher the bigger the c s parameter is, i.e. the less the 
concentration of the density of structure is. 
We, therefore, conclude that by interpreting astrophys- 
ical data of gravitational lensing by clusters using a fa- 
miliar NFW density profile without a void we are devi- 
ating from the exact result and the Th-L approximation 
is no longer valid. The Th-L approximation may, how- 
ever, be considered as precise enough if one modify the 
density profile and add the corresponding void to it, as 
require by general relativity for a quasi-local structure. 
The detail of the void, such as its density contrast, its 
depth and length, depends on the detail of the model 
and the deviation from the familiar NFW may even be 
much higher for other choices. Also note that the effect 
of the void is higher for larger impact parameter. In the 
case of strongly lensed objects in astrophysical applica- 
tions we are usually faced with small impact parameter 
where this effect is negligible. For example in the case of 
Abell 2261 cluster (z — 0.225) with many strong lensing 
arcs, D. Coe et al. [27! have assigned c s = 6.2 ± 0.3 and 
M mr = 2.2±0.2 x 10 15 M o . The exact general relativistic 
results according to our model would lead to c s = 6.23 
and M vir — 2.23 x 10 15 M Q . In the case of weak lensing, 
however, we expect this effect to have significant impact 
on the cosmological parameters. Work in this direction 
is in progress. 
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FIG. 6: Deviation angle for NFW density profiles with dif- 
ferent parameters. Horizontal axis is normalized to R s and 
vertical axis is normalized to the maximum of the deflection 
angle in each case. Dash line is for NFW model without void 
(formula |l9|)). 
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